clear all
set more off, perm
local mydir "`1'"
global Tables 	`mydir'/TS/Tables
global Work   	`mydir'/TS/DataWork
global Raw    	`mydir'/TS/DataRaw
global WRDS  	`mydir'/WRDS
/* ********************************************************************* */
import excel $Raw/shiller_updated.xls, sheet("formatted") firstrow clear
keep if Date<.
drop Fraction

foreach var of varlist RealDiv* RealEar* CAPE CAPD RCAE RCAD D E {
	destring `var', replace
}
* generate Stata dates
gen yr=int(Date)
gen mo=100*(Date-yr)

generate statpers=mdy(mo+1,1,yr  )-1 if mo<=11
replace  statpers=mdy(1   ,1,yr+1)-1	if mo==12
format statpers %td 

gen datem=mofd(statpers)
format datem %tm
* ******************************
* ******************************
keep if statpers<=mdy(2,1,2021)
* ******************************
* ******************************
keep datem RealEarnings RealDividend RealPrice RCAE CPI CAPE CAPD D E P
keep if CPI<.
sort datem

gen deflator=CPI/CPI[_N]

tempfile shiller
save "`shiller'", replace
* ******************************************************************************
use datem statpers fcast* LTG LTG_Nagel spindx if LTG<. using $Work/prgm02, clear
keep if statpers<=mdy(1,1,2021)
replace LTG=LTG/100
tsset datem

sort datem
ipolate fcastE1 datem, gen(ibesE)
ipolate fcastE5 datem, gen(ibesE5)
/* *************************** */
* add Shiller data for the period that overlaps with IBES --> * r and g are computed using data IBES period
merge 1:1 datem using "`shiller'", keepusing(datem deflator D E P) update
keep if _merge==3
drop _merge
/* *************************** */
* note: r and g are computed using data IBES period
scalar r=(spindx[_N]/spindx[1])^(12/_N)-1
scalar rm=(1+r)^(1/12)-1

scalar g=(D[_N]/D[1])^(12/_N)-1

gen DE=D/E
quietly: su DE
scalar pay=r(mean)

gen de = ln(D/E)
quietly: su  de 
scalar m_de=r(mean)

gen dp = ln(D/P)
quietly: su  dp 
scalar rho=1/(1+exp(r(mean)))
drop DE de dp 

merge 1:1 datem using "`shiller'"
drop _merge
* ******************************************************************************
* add Shiller data for observations outside of the IBES range
merge 1:1 datem using "`shiller'", update 
drop _merge

tsset datem
gen STG1_E=ln(fcastE1/L3.E)
gen STG2_E=ln(fcastE2/fcastE1)

gen STG1_D=ln(fcastD1/L3.D)
gen STG2_D=ln(fcastD2/fcastD1)

* Shiller's data on dividends and earnings ends on 12/2019
keep if E<.
* ******************************************************************************
gen p=ln(P)
gen d=ln(D/12)
gen e=ln(E/12)

scalar alfa=rho
scalar k  = -ln(alfa ) - (1-alfa ) * ln(1/alfa -1)
scalar ke = k + (1-alfa) * m_de

scalar alfam=rho^(1/12)  				/* page 134 Campbell's book. bottom paragraph	*/
scalar km = -ln(alfam) - (1-alfam) * ln(1/alfam-1)	/* page 134 Campbell's book. equation 5.40	*/
scalar ke = km + (1-alfam) * m_de			/* De la O and Myers (page 9)	 		*/

* new pRE index
gsort -datem

generate pRE_D=(km-rm)+alfam*ln(1.0*D[_n-1]/(r-g))+(1-alfam)*d[_n-1] 		if _n==2
generate pRE_E=pRE_D 								if _n==2

replace  pRE_D=(km-rm)+alfam*pRE_D[_n-1]+(1-alfam)*d[_n-1] 			if _n>2
replace  pRE_E=(ke-rm)+alfam*pRE_E[_n-1]+(1-alfam)*e[_n-1] 			if _n>2
/* *************************** */
tsset datem
* compute price index using dividends growth
gen pv1=alfa^0*ln(fcastD1/D)+ln(D)+(ke-r)/(1-alfa)
gen pv2=alfa^1*ln(fcastD2/fcastD1)+pv1

gen ni1=pv1
gen ni2=pv2

forvalues x=3(1)50 {
	local j=`x'-1
	gen pv`x'=alfa^`j'*ln(1+LTG)+pv`j'

	local w=int((`x'-1)/5)
	gen ni`x'=alfa^`j'*ln(1+0.4^`w'*LTG)+ni`j'
}
*************************
gen pv1_v=alfa^0*ln(fcastD1/D)+ln(D)+(ke-r)/(1-alfa)

forvalues x=2(1)50 {
	local j=`x'-1
	gen pv`x'_v=alfa^`j'*ln(1+LTG)+pv`j'_v
}
forvalues x=1(1)50 {
* terminal growth rate fits average prices
	gen g`x'=(1-alfa)*(p-pv`x'_v)/alfa^`x'
	egen gm`x'=mean(g`x')
	gen p`x'_D_v=pv`x'_v+gm`x'*alfa^`x'/(1-alfa)

* sanity check (make sure that i am recovering the actual stock price)
	gen m`x'=pv`x'_v+g`x'*alfa^`x'/(1-alfa)
}
order *, sequential
drop m1-m50 g1-g50 gm1-gm50 
*************************
gen Adi=pv10+k/(1-alfa)-p

tsset datem
quietly: su LTG
scalar mLTG=r(mean)

forvalues x=1(1)50 {
* terminal growth rate fits average prices
	gen g`x'=(1-alfa)*(p-pv`x')/alfa^`x'
	egen gm`x'=mean(g`x')
	gen p`x'_D=pv`x'+gm`x'*alfa^`x'/(1-alfa)

* sanity check (make sure that i am recovering the actual stock price)
	gen m`x'=pv`x'+g`x'*alfa^`x'/(1-alfa)

* terminal growth rate fits average prices (nicola's experiment)
	drop g`x' gm`x'
	gen g`x'=(1-alfa)*(p-ni`x')/alfa^`x'
	egen gm`x'=mean(g`x')
	gen nD`x'=ni`x'+gm`x'*alfa^`x'/(1-alfa)
}

pwcorr m1 m2 m50 p

order *, sequential
drop m1-m50 g1-g50 gm1-gm50 pv* ni*
* **************************************
* work with earnings
scalar ke= k + (1-rho) * ln(pay)		/* De La O and Myers pag 9 */
gen pv1=alfa^0*ln(fcastE1/E)+ln(E)+(ke-r)/(1-alfa)
gen pv2=alfa^1*ln(fcastE2/fcastE1)+pv1
gen pv3=alfa^2*ln(fcastE3/fcastE2)+pv2
gen pv4=alfa^3*ln(fcastE4/fcastE3)+pv3
gen pv5=alfa^4*ln(fcastE5/fcastE4)+pv4

gen ni1=pv1
gen ni2=pv2
gen ni3=pv3
gen ni4=pv4
gen ni5=pv5

forvalues x=6(1)50 {
	local j=`x'-1
	gen pv`x'=alfa^`j'*ln(1+LTG)+pv`j'

	local w=int((`x'-1)/5)
	gen ni`x'=alfa^`j'*ln(1+0.4^`w'*LTG)+ni`j'
}

forvalues x=1(1)50 {
* terminal growth rate fits average prices
	gen g`x'=(1-alfa)*(p-pv`x')/alfa^`x'
	egen gm`x'=mean(g`x')
	gen p`x'_E=pv`x'+gm`x'*alfa^`x'/(1-alfa)

* sanity check (make sure that i am recovering the actual stock price)
	gen m`x'=pv`x'+g`x'*alfa^`x'/(1-alfa)

* terminal growth rate fits average prices (nicola's experiment)
	drop g`x' gm`x'
	gen g`x'=(1-alfa)*(p-ni`x')/alfa^`x'
	egen gm`x'=mean(g`x')
	gen nE`x'=ni`x'+gm`x'*alfa^`x'/(1-alfa)
}
order *, sequential
pwcorr m1 m50 p

drop m1-m50 g1-g50 gm1-gm50 ni*

rename p lnp

order *, sequential

tsset datem
foreach var of varlist lnp pRE_D p1_D p2_D p10_D nD10 nD15 nD20 {
	generate d`var'=`var'-L12.`var'
	generate `var'd=`var'-ln(D)
}	
drop dlnp
foreach var of varlist lnp pRE_E p1_E p2_E p10_E nE10 nE15 nE20 {
	generate d`var'=`var'-L12.`var'
	generate `var'e=`var'-ln(E)
}	
drop if dofm(datem)>=mdy(1,1,2021)
* ****************************
save $Work/figure01, replace
* ****************************